An Investigation of Mean-field Effects for a Bose Condensate in an Optical Lattice 



S. B. McKagan 1 '*, D. L. Feder 2 , and W. P. Reinhardt 1 
1 Departments of Chemistry and Physics, University of Washington, Seattle, WA 98195 and 
2 Department of Physics and Astronomy and the Institute for Quantum Information Science, 
University of Calgary, Calgary, Alberta, Canada T2N 1N4 
(Dated: February 2, 2008) 

This paper presents a mean-field numerical analysis, using the full three-dimensional time- 
dependent Gross-Pitaevskii equation (GPE), of an experiment carried out by Orzel et al. [Sci- 
ence 291, 2386 (2001)] intended to show number squeezing in a gaseous Bose-Einstein condensate 
in an optical lattice. The motivation for the present work is to elucidate the role of mean-field 
effects in understanding the experimental results of this work and those of related experiments. We 
show that the non-adiabatic loading of atoms into optical lattices reproduces many of the main 
results of the Orzel et al. experiment, including both loss of interference patterns as laser intensity 
is increased and their regeneration when intensities are lowered. The non-adiabaticity found in the 
GPE simulations manifests itself primarily in a coupling between the transverse and longitudinal 
dynamics, indicating that one-dimensional approximations are inadequate to model the experiment. 

PACS numbers: 03.75.Dg, 03.75.Lm, 05.30.Jp, 32.80.Pj 



I. INTRODUCTION 

The creation of dilute gaseous Bose-Einstein conden- 
sates (BECs) in the laboratory in 1995 [1 S H has 
spurred much development in both experiment and the- 
ory 0. 0. Q Mean- field theory was able to explain most 
early BEC experiments, using the well-known Gross- 
Pitaevskii equation (GPE) 0,15 which is valid for weakly 
interacting BECs at zero temperature because it assumes 
all atoms occupy a single macroscopic wavefunction: 

^ft^(r, t) = (-t^V 2 + V(r, t) + g^r, i)| 2 ) V(r, t) 

(1) 

where tp(r,t) is the single particle wave function for any 
atom in the BEC, V(r,t) is an external trapping poten- 
tial, and g = Aita, s h 2 /m, where a s is the s-wave scattering 
length and m is the atomic mass. These early successes 
included the anisotropic profile and momentum distribu- 
tion of the ground state after ballistic expansion 0, U ' 
the spectrum of collective excitations |a |9j , the dynamics 
of spinor condensates 0, Ell . and the wave interference 
of interacting condensates |l2l E^ |. Perhaps unexpect- 
edly, the time dependent GPE also accurately described 
strongly nonlinear excitations such as vortices p^ . EolE^| 
and solitons [TtLEH . 

Given the overwhelming success of the GPE, there 
has been much interest in finding situations in which 
it breaks down and a more detailed theoretical descrip- 
tion is needed. Indeed, the theoretical interpretation 
of the earliest experiments on collective excitations at 
finite temperature |9| Il9l |20j has required a dynami- 
cal theor y t hat includes the motion of the nonconden- 
sate 0, [2^ ■ More recent experiments are geared to- 
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ward inducing strong correlations among atoms, which 
are not captured by the GPE- these include molecules 
in Bose [23, 01 an d Fermi I25l.l26l.l27l gases and explor- 
ing the BCS-BEC crossover |28|,|29, 30], low-dimensional 
systems through tight confinement |3ll 133.1331 . quantum 
Hall-like states in rapidly rotating traps |34Ll35j or using 
external lasers [36 | . an d novel many-body states of atoms 
in optical lattices |Hl |H |H S E3, El ■ 

In the experiments with optical lattices, the extent of 
correlations in the atomic gas has generally been mea- 
sured by dropping all confining potentials and imaging 
the cloud after a period of ballistic expansion. When 
condensates are suddenly released from their confine- 
ment in shallow optical lattices formed from weak lasers, 
they form well-defined interference patterns correspond- 
ing to momentum-space Bragg peaks. This has been 
interpreted as an unequivocal signature of phase coher- 
ence over multiple lattice sites. As the laser intensity 
increases, however, the interference patterns partially 
or fully wash out [37], I2E H^]. The theoretical inter- 
pretation is that as tunneling is quenched, the initial 
macroscopic condensate splits into many separate sub- 
condensates which may then lose some or all of their rel- 
ative coherence; the resulting state has been described as 
'number squeezed,' or 'fragmented.' Under various con- 
ditions one obtains full fragmentation, which leads to a 
quantum phase transition from a superfluid state to a 
Mott insulator 0. 

In fact, the loss of fringe contrast after ballistic expan- 
sion is not as clear a signature of condensate fragmen- 
tation as is often assumed. Jean Dalibard and collabo- 
rators |44| recently demonstrated that high-visibility in- 
terference patterns resulted even under conditions where 
the phases from site to site of a one-dimensional (ID) 
optical lattice were random. A simple theoretical model 
that generalizes the interference pattern from two uncor- 
related sources reproduces the experimental data. Under 
conditions similar to those of Orzel et al. 381 where num- 
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ber squeezing is expected to be present (though using 
a blue-detuned rather than a red-detuned lattice where 
the transverse confinement is much weaker), they ob- 
served an unattributed heating effect. A similar heat- 
ing was observed by Morsch et al. 45], where radial 
modes excited by the non-adiabatic ramp-up of the ID 
optical lattice rapidly damped and transfered energy to 
high- lying axial modes. Zakrzewski |46j|. using a time- 
dependent Gutzwiller mean-field approach to solving the 
Bose-Hubbard model, has found that non-adiabatic mean 
field effects reproduce some of the results seen in the ex- 
periments of Greiner et al. |39j |. 

The recent experiments indicate that fundamental 
questions need to be answered before fringe contrast can 
be used as a quantitative measure of the breakdown of 
mean-field theory in these systems. The effects of non- 
adiabaticity on the phase variations across the lattice 
have been addressed previously 0, 0> US, ■ ID cal- 
culations based on the GPE found that the harmonically 
trapped condensate develops a pronounced 'phase sag' af- 
ter slowly ramping up the lattice 0, 0] ; the phase was 
found to increase approximately quadratically around 
the trap center. A flat phase profile at the end of the 
ramp, and a high-contrast interference pattern, can be 
restored if the harmonic trap parameters are simultane- 
ously varied |4^. More recent effective 3D calculations 
demonstrate that collective modes are excited by slow 
ramps 0, OH] , but the influence of these excitations on 
the resulting interference patterns was not addressed. 

In order to clarify the various mean-field effects that 
can degrade the interference pattern when a condensate 
is released from an optical lattice, we have performed 
fully 3D numerical simulations based on the GPE of one 
of the above experiments, namely that of Orzel, Tuch- 
man, Fenselau, Yasuda, and Kasevich (33 (referred to 
as OTFYK in what follows). See also |52|. If number 
squeezing were the main reason for the observed loss of 
fringe contrast as has been suggested |3^| , the GPE sim- 
ulations should not be able to mimic the experimental 
data. We find, however, that we can reproduce the degra- 
dation of the interference patterns as the lattice depth ap- 
proaches the regime where squeezing would be expected, 
and the subsequent restoration of contrast as the laser 
intensity is lowered, even in cases where 'random' phases 
have been applied to the individual wells at the largest 
lattice depths. At the level of the mean-field approxima- 
tion, the behavior is due to the non-adiabatic response 
of the trapped condensate as the optical lattice is slowly 
turned on: the ramp induces axial currents and the ex- 
citation of strongly coupled transverse and longitudinal 
excitations. These in turn yield variations in the phase 
from site to site that cause the interference patterns to 
disappear. We expect that in the presence of damping 
(not included in the present analysis) , the mean- field re- 
sults would be close to, if not indistinguishable from, the 
experimental data. 

The paper is organized as follows: In Section II, the 
ground states are obtained for a condensate in harmonic 



oscillator, gravitational, and optical lattice potentials as 
a function of laser intensity. The OTFYK experiment is 
described in Section III, and the observations are com- 
pared to results of a GPE model that assumes perfect 
adiabaticity of the optical lattice loading. Section IV, 
which is the heart of the paper, contains discussion of 
a GPE simulation of the full experimental protocol, in- 
cluding the turning on and turning off of the lattice, and 
the use of gravity to induce a relative phase shift of 7r 
between adjacent lattice sites. It is found that the ex- 
perimental timescales for the lattice ramp are not gener- 
ally sufficient to ensure adiabaticity. The resulting phase 
shifts mimic loss of coherence between adjacent lattice 
sites. In Sections V and VI it is demonstrated that in- 
terference patterns recover as the laser intensities are re- 
duced, even if the site-to-site phases are artificially ran- 
domized at lattice maximum. In addition to the accumu- 
lation of axial phase variations, a second source of non- 
adiabaticity is found and discussed in Sections VII and 
VIII: a strong coupling between longitudinal and trans- 
verse oscillations. The results are compared to more re- 
cent experiments by the Kasevich group [U 01 • A brief 
summary and conclusions end the paper, in Section IX. 

II. THE POTENTIAL AND GROUND STATE 
DENSITIES 

In the OTFYK experiment, a BEC is first created in 
a cylindrically symmetric magnetic trap with the longi- 
tudinal axis oriented vertically, parallel to the ambient 
gravitational field. Counter-propagating laser beams are 
then slowly turned on, oriented vertically and aligned 
with the long axis of the magnetic trap. The lasers are 
red detuned, and thus the condensate is pulled into the 
anti-nodes, with spacing A/2, and also experiences strong 
transverse confinement. The full potential used to de- 
scribe in the present numerical studies of this system is 
given by [H^: 

V{p, z, t) = U(t) (l - e - p2/r ?» cos 2 (fcz)) 

+ ^m{Lu 2 ±P 2 + u 2 z z 2 ) + mgz (2) 

where the first term is due to the laser, the second term 
to the magnetic trap, and the third term to gravity. The 
longitudinal coordinate z is in the vertical direction and 
the transverse coordinate p — y 1 x 1 + y 2 . Table [I] gives 
the values of the experimental parameters as reported in 
the OTFYK paper. 

Various parts of the trap are turned on and off at vari- 
ous points in the experiment, as described in Section III, 
so not all the terms in Eq. (J2J) are present at all times. 
The laser is generally ramped up slowly, so U varies as 
a function of time. As long as the magnetic trap is on, 
the gravitational term has no effect other than to give a 
linear shift the quadratic trapping potential, but it plays 
an important role when the magnetic trap is turned off. 
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TABLE I: Values of parameters used in the experiment and 
modeling. 



N 


30000 


number of atoms in BEC 


a s 


5.24nm [56] 


scattering length of 87 Rb 


m 


1.44 x 10 _25 kg [57] 


atomic mass of 87 Rb 


LU ± 


2tt x 19Hz 


transverse frequency of trap 






longitudinal frequency of trap 


k 


27r/840nm 


wave number of laser 


rib 


50/xm 


1/e intensity radius of laser 


U 


to 44E R 


Laser intensity (E R — fi 2 k 2 /2m) 



In our calculations, we neglect this term except when it 
is required for phase manipulation. Note that the laser 
creates not only a periodic trapping potential in the longi- 
tudinal direction, but also a Gaussian trapping potential 
in the transverse direction. This Gaussian term (which 
is not present for a blue-detuned laser) implies that as 
the laser strength increases, so does the transverse con- 
finement of the BEC, and therefore the density increases. 
This turns out to be a very important effect. 

Although the actual OTFYK experiments, as seen in 
Section III, consist of a sequence of well defined time de- 
pendent steps, it is useful to examine the density profiles 
of the stationary GPE ground states in the potential of 
Eq. Pjl. Ground states were obtained using numerical 
methods described in the Appendix, by evolution of the 
GPE in complex time (t — > it) with propagation converg- 
ing to the lowest energy stationary states for each fixed 
value of U. The results for three representative values of 
U are shown in Fig-QJ Figure^a) illustrates the density 
profile in the absence of the optical lattice. All experi- 
ments and theoretical simulations described herein start 
with this ground state of the condensate in the magnetic 
trap in the presence of gravity. Figures ^b) and (c) il- 
lustrate the two dominant effects of the laser fields: first 
the condensate morphs from a single ellipsoidal and co- 
herent structure into a vertical stack of disk like sub- 
condensates. Also clearly seen is the very large effect of 
the transverse confinement: the condensate density en- 
velope changes from oblate to prolate at the highest field 
intensities. 



III. AN INITIAL "ADIABATIC" SIMULATION 
OF THE EXPERIMENT 

We have simulated each step of the OTFYK experi- 
ment using the full 3D GPE with the potential discussed 
in the previous section. To avoid repetition, we will de- 
scribe the steps of the experiment and our simulation in 
parallel. The first step of the experiment is to create a 
BEC in a magnetic trap. This step can be simulated by 
using complex time evolution as described in the previous 
section. The results of the simulation for the harmonic 
trap only are shown in Fig. flfa) . 

The next step in the experiment involved turning on 
the laser fields so that the intensity increases linearly 
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FIG. 1: Cross sections of integrated density for the ground 
state of a BEC in a harmonic trap and laser, from complex 
time evolution. Note that the effect of the transverse confine- 
ment of the laser is very large, turning an initially disk-shaped 
BEC into a cigar-shaped BEC. 



from zero to some final intensity Uf (ranging from 7.2Er 
to AAEr) in a ramp time tr = 200 ms. This ramp- 
up time was assumed by OTFYK to be slow enough to 
allow the condensate to follow the ramping adiabatically, 
so that the BEC stays in the appropriate ground state 
as the laser is turned up [H^; in the next section we will 
test this assumption of adiabaticity. If the final state 
of the BEC in the combined magnetic trap and lattice 
is the ground state of this system, then the sequence of 
states generated are just those which may be found by 
complex time evolution, as discussed in Section II. Then 
the ground states as a function of U shown in Figs.^b,c) 
can be used in a preliminary simulation of the OTFYK 
experiment. 

The OTFYK paper described two regimes that exhib- 
ited markedly different interference patterns after the ex- 
ternal potentials were dropped and the atomic cloud was 
allowed to ballistically expand. For low Uf, the final 
state was expected to be fully phase coherent yielding a 
clean interference pattern jH^] . The computed three-peak 
pattern for the Uf = 7.2Er case is shown in Fig. |5Ja). 
For high Uf ~ AOEr, the condensate was anticipated to 
be strongly number squeezed, leading to a random rel- 
ative phase from site to site. Ballistic expansion would 
then yield no detectable interference pattern. Of course, 
the GPE cannot produce a number squeezed state, but 
if the random phases are put in 'by hand' it can mimic 
the experimental situation, as was done previously for 
the MIT interference experiment 01 by Rohrl et al. 0] 
within the framework of the GPE. The numerical results 
for this case are shown in Fig. H[b). 

The three peak interference pattern of Fig. HJa) is in- 
convenient for quantitative estimation of the extent of 
decay of coherence, so a third protocol was invoked by 
OTFYK. The experimenters used the presence of gravity, 
via an appropriate time delay, as discussed in the follow- 
ing section, to imprint a 7r phase difference between the 
condensate in each pair of neighboring wells. This pro- 
duces a two-peak interference pattern, as illustrated in 
the simulated ballistic interference pattern of Fig. |2Ic) . 
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FIG. 2: Numerical simulations of three peak, washed out, 
and two peak interference patterns produced by imprinting 
phases on ground state of GP wave function in harmonic trap 
and 7.2Er laser. The first row shows the phase along the 
longitudinal axis before expansion. The second row shows 
the density along the longitudinal axis after releasing the BEC 
and allowing it to expand for 8 ms. The third row shows a 
cross section of the density after expansion (indicated by the 
gray scale). Note the difference in the scale before and after 
expansions. 



The three situations of Fig. [21 were simulated by paint- 
ing phase patterns 'by hand' onto the ground state GP 
wave function in a 7.2Er lattice and then allowing it 
to spatially expand by turning off all external potentials. 
The next section will show the condensate phase patterns 
obtained by a full time-dependent simulation of the ac- 
tual experiment. 



IV. FULL TIME-DEPENDENT GPE 
SIMULATION OF THE OTFYK EXPERIMENT 

We have tested the assumption of adiabaticity by simu- 
lating a 200 ms laser turn-on in real time. In all cases the 
optical lattice was turned on starting with a numerically 
exact ground state condensate in the magnetic trap. The 
results of these tests, summarized in Fig. |3 are that a 
200 ms switching time is nearly adiabatic for Uf = 7.2Er, 
but not for Uf = 40Er. Much longer times are needed 
for higher laser intensities not only because the laser is 
turned to higher intensity, but because the higher laser 
fields push more of the condensate into the outer wells; 
the required times for the condensate to tunnel through 
higher barriers to reach these outer wells are also longer. 

As the lattice is raised, the value of the chemical po- 
tential varies from site to site, and atoms must flow 
from the trap center to the periphery in order to remain 
in equilibrium. For sufficiently deep lattices, the time 
needed to tunnel from site to site eventually exceeds the 
timescale of the ramp. The resulting non-adiabaticity 
for the Uf — AOEr case leads to significant deviations of 
both the density and phase profiles, compared with the 
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FIG. 3: Density profiles of ground state with the trap and 
lattice compared with density profiles after turning up laser. 
Each plot shows the density and phase along the longitudinal 
axis and z > 0. The first row shows the density profiles of 
the ground state of the condensate in the trap and lattice for 
U — 7.2Er and 40Er. The second row shows the density 
profile after creating the condensate in the trap and turning 
up the laser in tr — 200 ms for the same laser strengths. If the 
laser turn-on is adiabatic, the profiles in the first and second 
rows should be identical. The profiles are nearly identical for 
7.2Er, aside from a very small phase sag, showing that the 
turn-on is adiabatic for this laser strength. For 40Er after a 
200 ms turn-on, however, the density profile is distorted at the 
edges and the phase sag cycles over more than 6n, showing 
that a longer turn on time would be required for adiabaticity 
at this laser strength. The third row shows the density profile 
after turning the laser up to 40Er in 2 s. In this case, the 
density matches that of the ground state and the phase sag is 
much smaller, showing that this turn-on is nearly adiabatic. 



true ground state. First, the density envelope is trun- 
cated, reflecting the inability of atoms to fully tunnel 
out to the cloud surface. Second, there is considerable 
'phase sag,' illustrated in the second column of Fig. 
This phase profile, which results from the axial veloci- 
ties acquired by the atoms as they propagate (v oc V<p), 
eventually becomes locked in for very deep lattices. 

A simple estimate of the timescale required for adi- 
abaticity is that it should at least exceed the inverse 
of the smallest collective mode frequency in the pres- 
ence of the lattice. In deep lattices where the atomic 
profile in a given well is not much different from that 
of an ideal gas, the effective axial frequency shifts as 
Cu z = lu z \Jmjra* , where m* is the effective mass of the 
atom 60]. For a 40E R lattice, m*/m w 900 HJ, which 
yields 2it/uj z w 600 ms. It is reasonable to expect the full 
adiabatic ramp to require several times this, tn ~ 2 s, 
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as is fully confirmed by the numerical results shown in 
Fig. |21 This value is an order of magnitude larger than 
that used in the experiment; since it is comparable to the 
lifetime of the BEC, such a long ramp is probably not 
experimentally feasible. An alternative protocol would 
have been to grow a coherent condensate ground state in 
the presence of the magnetic trap and optical lattice. 

After turning the laser up to Uf in 200 ms, the experi- 
menters turn off the magnetic trap in 40 /its and hold the 
BEC in the vertically oriented laser for approximately 
2.5 ms in order to produce a it offset between adjacent 
wells. The confinement of the laser is sufficient to prevent 
significant movement of the condensate during this time. 
However, the gradient of the gravitational potential be- 
tween neighboring wells of the lattice causes dramatic 
Schrodinger phase evolution, since the condensate phase 
is given by — Vt/h. Since there is a difference of po- 
tential energy between the wells equal to AV^ = mgAz, 
where Az — A/2 is the distance between the wells, after 
the condensate is held in the laser for a time th there 
will be a phase difference between the condensates in 
two neighboring wells equal to A8 = mg\th/2h. To 
achieve AO = 7r, the hold time must be an odd multi- 
ple of 2ith/mg\ — 0.557ms for the parameters used in 
this experiment. 

In our calculations, we follow this experimental pro- 
cedure, using the exact hold time for AO = 5ir, th = 
2.785 ms, rather than the 2.5 ms quoted in OTFYK |6lj . 
Simulations with th = 2.5 ms give an asymmetric density 
distribution after expansion, in which one peak is about 
twice as large as the other. This asymmetry is due to 
the asymmetry of the Fourier components of the wave 
function when the phase difference between the wells is 
not exactly it [37|. 

A. Real Time Evolution for a Weak Laser 



After tumina on laser in 200 ms? 
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FIG. 4: Real time evolution for a weak laser, (a) Density 
profile along the longitudinal axis after turning laser up to 
7.2Er in 200 ms. (b) Density profile along the longitudinal 
axis after turning off the magnetic trap holding the BEC in 
the laser (with gravity) for 2.785 ms. (c) Phase profile along 
the longitudinal axis after turning laser up to 7.2Er in 200 ms. 
(d) Phase profile along the longitudinal axis after turning off 
the magnetic trap holding the BEC in the laser (with gravity) 
for 2.785 ms. (e) Density profile along the longitudinal axis 
after releasing the BEC and allowing it to expand for 8 ms. 



For a weak laser, the 200 ms turn-on is nearly adia- 
batic, but it is still important to check the results of the 
phase evolution and expansion of the actual state after 
the laser is turned on. The results of this calculation 
are shown in Fig. These results are as expected in 
that there is approximately a it phase difference between 
each neighboring well, and a clear two-peak interference 
pattern after expansion. Thus, although the simulated 
dynamics of the BEC in the weak laser are not as clean 
as the predictions of Fig. [3 the simulations confirm that 
the basic ideas of the discussion of Section III are correct 
in this case, and in agreement with the OTFYK obser- 
vations. 



B. Real Time Evolution for a Strong Laser 

Fig.[5]illustrates the effect of ramping the lattice up to 
40.Br in 200 ms. The resulting phase sag shown in the 
second column of Fig. [21 is repeated in Fig. [SJc). After 



holding the atoms in the laser for 2.785 ms, as shown in 
Fig. |^d) , the phase is so distorted that it begins to re- 
semble the random pattern used in Fig. [2 In this context 
it is not surprising that the density profile after expan- 
sion, depicted in Fig.EIe), is completely washed out. The 
loss of interference is driven entirely by non-adiabatic ef- 
fect captured within a mean-field model. It is important 
to underline that these non-adiabatic effects are intrinsic 
to the ramp time, and are not affected by how the ramp 
is applied. We have performed simulations using a ramp 
with a smooth onset (based on a sine function) rather 
than the linear ramp discussed above, but the observed 
loss of interference was unchanged. We are thus able to 
qualitatively reproduce the loss of interference without 
invoking number squeezing. 
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FIG. 5: Real time evolution for a strong laser, (a) Density 
profile along the longitudinal axis after turning laser up to 
40Er in 200 ms. (b) Density profile along the longitudinal 
axis after turning off the magnetic trap holding the BEC in 
the laser (with gravity) for 2.785 ms. (c) Phase profile along 
the longitudinal axis after turning laser up to 40Er in 200 ms. 
(d) Phase profile along the longitudinal axis after turning off 
the magnetic trap holding the BEC in the laser (with gravity) 
for 2.785 ms. (e) Density profile along the longitudinal axis 
after releasing the BEC and allowing it to expand for 8 ms. 

C. Absorption Probabilities 

To compare the results of our simulations directly with 
the experimental data, we calculated absorption proba- 
bilities, which we then smoothed to account for finite 
experimental resolution. For each simulation, we inte- 
grate the 3D density profile after 8 ms of expansion, 
N\ip(x, y, z)\ 2 , over one of the transverse axes (the direc- 
tion of the imaging beam), to get the integrated density: 

(3) 

We then calculate the absorption probability from the 
following equation 62]: 

A{x, z) = 1 — exp (— crofi(x, z)) (4) 

where the absorption cross section tro is given by: 

(Tq = 67rA 2 (5) 



with A = 780nm/27r. We then smooth the data by taking 
a Gaussian convolution: 

A(x, z) = J dx'dz'A{x', z >) e -((*-*') 2 +(z-z') 2 )/™ 2 ( 6) 

where the 1/e width w — 17//m is chosen so that at very 
low laser strength (Uf — 5 — 8 En) the ratio of the width 
of the interference peaks to their separation is £ = 0.22. 
This value of £ corresponds to the best experimental con- 
trast reported in Ref. [38j (the parameter £ will be dis- 
cussed further in Sec. lVIIl|l . The same method was used 
by OTFYK to analyze their experimental results, but 
they found that they needed a 1/e width of 25/im, rather 
than 17/Km, to match their data |63j . 

Fig.Elshows how the results of time-dependent numer- 
ical simulations compare with the experimental results. 
The three columns depict the raw numerical results for 
integrated density, the smoothed absorption probabili- 
ties, and the experimental results (copied with permis- 
sion from the authors), respectively. 

The numerical simulations can reproduce qualitatively 
the observed interference patterns for both weak and 
strong lasers. The fiat phase profile after loading the 
atoms into a 7 ' .IEr lattice guarantees a clean interfer- 
ence pattern |59| . For a strong laser, the underlying 
mechanism for the loss of interference in our simulations 
(phase distortion due to non-adiabatic mean field effects) 
is entirely different from the mechanism proposed by the 
experimenters (number squeezing). A more quantitative 
analysis, given in Section VIII, reveals that the corre- 
spondence between simulation and experiment is not per- 
fect, so it may be that number squeezing is in fact oc- 
curring in the experiment. However, the results of these 
simulations show that loss of interference is not in itself 
sufficient evidence of the presence of number squeezing. 



V. RESTORING INTERFERENCE 

In the next phase of the experiment, instead of releas- 
ing the BEC after turning the laser up to 40_Er, OT- 
FYK first turn the laser down to IOEr in 150 ms and 
then release the condensate after the short hold time. 
In the experiments the interference pattern, and by im- 
plication the coherence of the BEC, is restored by the 
adiabatic ramp-down. The results were originally ex- 
plained in terms of the time-dependent two-mode model 
as the loss and return of coherence of the BEC. How- 
ever, Fig. [7| illustrates that they can also be reproduced 
with numerical simulations of the GPE, and can there- 
fore be explained in terms of mean-field effects. Indeed, 
the original interpretation of these experiments in terms 
of number squeezing has since been revised (5j| . The dy- 
namics are illustrated in Fig. El which shows the density 
and phase profiles after turning down the laser, before 
and after releasing the BEC. This figure shows that af- 
ter turning the laser up and then down, the phase sag 



7 



Simulations 

Smoothed 
Integrated Absorption 
Densities Probabilities 




Experiments 

Absorption 
Images 



O 



Simulated Integrated Densities! 



B 



C 



t 



» 



50 



FIG. 6: Comparison between experimental results and nu- 
merical simulations. Interference patterns after turning laser 
intensity up to Uf in 200 ms, turning off trap and holding 
BEC in laser for 2.785 ms, then releasing it and allowing it 
to expand for 8 ms. (a) U f = 7.2E R (b) U f = 18E R (c) 
u f = AAE R 
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FIG. 7: Restoring interference: Comparison between numer- 
ical and experimental results, (b) Interference pattern after 
releasing BEC from 7.2E R laser, (c) Interference pattern after 
turning laser up to A0E R and releasing BEC. (d) Interference 
pattern after turning laser up to AOEr, then turning laser 
down to 10-Er, and then releasing BEC. 



unwinds, so that at the end of the process, the site-to- 
site phase is relatively constant. There are variations in 
phase, but their size is on the order of 7T, rather than 77r. 
There is considerably more noise in the system after this 
process than there would have been if the laser were sim- 
ply turned up to IOEr without first going through the 
high barrier state, but this noise does not obscure the 
basic pattern. 



VI. RANDOM PHASE IMPRINTING 

A further step, which was not done in the experiment, 
but which in principle could be done, is to imprint a ran- 
dom phase shift on each well when the laser strength is 
AOEr and then turn the barrier back down. The two- 
mode model predicts that applying a random phase shift 



will destroy the ability of a mean-field state to heal back 
to a stationary (flat-phase) BEC state, but will have no 
effect on a strongly number-squeezed state. Since we 
have shown that the loss of interference is not sufficient 
to demonstrate number squeezing, it appears that ran- 
dom phase shifts could be used as an alternative test. If 
this is to be an effective test, applying a random phase 
shift should completely destroy the ability of the result- 
ing BEC to produce a clean interference pattern. 

We tested this idea numerically by running simulations 
similar to those described in the previous section, where 
we turned the laser up to AOEr and then down to WEr, 
but in this case we applied a random phase shift to each 
well before turning the laser down. The surprising result, 
illustrated in Fig. [51 is that the non-linear mean field dy- 
namics 'self-heals' a truly random phase pattern just as it 
unwinds a mean-field induced phase sag: the interference 



8 



After turning laser down in 1?0 ins 
(a) (b) 




-8-4048 -8-4048 
z (|tira) z (|tini) 



After expanding for 8 ins 




-100 -50 50 100 

z (jj,m) 



FIG. 8: Real time evolution after turning laser up and then 
down, (a) Density profile along the longitudinal axis after 
turning laser up to 40Er in 200 ms and then down to IOEr 
in 150 ms. (b) Density profile along the longitudinal axis 
after turning off the magnetic trap holding the BEC in the 
laser (with gravity) for 2.785 ms. (c) Phase profile along the 
longitudinal axis after turning laser up to 40Er in 200 ms and 
then down to IOEr in 150 ms. (d) Phase profile along the 
longitudinal axis after turning off the magnetic trap holding 
the BEC in the laser (with gravity) for 2.785 ms. (e) Density 
profile along the longitudinal axis after releasing the BEC and 
allowing it to expand for 8 ms. 



After turning laser down in 1?0 ms 
(a) (b) 
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FIG. 9: Real time evolution after turning laser up, applying a 
random phase shift to each well, and then turning laser down, 
(a) Density profile along the longitudinal axis after laser is 
turned down, (b) Density profile along the longitudinal axis 
after turning off the magnetic trap holding the BEC in the 
laser (with gravity) for 2.785 ms. (c) Phase profile along the 
longitudinal axis after turning laser up to 40Er in 200 ms and 
then down to 10-Ei? in 150 ms. (d) Phase profile along the 
longitudinal axis after turning off the magnetic trap holding 
the BEC in the laser (with gravity) for 2.785 ms. (e) Density 
profile along the longitudinal axis after releasing the BEC and 
allowing it to expand for 8 ms. 



VII. PHASE OSCILLATIONS 



pattern is nearly as dramatic as in the simulations in the 
previous section where there was no phase shift. In Fig.El 
the two-peak interference pattern is not very clean, but 
Fig.llOlshows that after smoothing, the absorption prob- 
ability is not significantly different from that acquired af- 
ter turning the laser up and down without a phase shift. 
Thus, even in the case of random phase shifts, mean field 
effects are able to mimic the predictions of the two-mode 
model and it appears to be difficult, under the OTFYK 
conditions, to distinguish experimentally between mean 
field effects and number squeezing. 



An unexpected effect observed in the numerical simu- 
lations during the turn-on of the laser is a slow oscillation 
of the phase. The phase sag does not continue to grow, 
as has been seen in previous simulations using effective 
ID models ^3 C3 > an d which one might expect if the 
phase dynamics were due only to the mean field poten- 
tial differences between the wells. Instead, as shown in 
Fig. ^] the phase oscillates, with the sag growing, then 
shrinking, then reversing. 

These phase oscillations are the result of two non- 
adiabatic effects. First, if the lattice is ramped up too 
quickly, the local chemical potentials /ii oc in each well, 
defined by the expectation of the GPE operator, will 
not be identical. At the end of the ramp, the phase in 
each well will vary in time as <^>i oc ~ H\ oc t/K, causing 
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FIG. 10: Integrated density and absorption probability of ex- 
panded BEC after turning laser up, applying a random phase 
shift to each well, and then turning laser down 

U(E R ) 
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FIG. 11: Phase accumulation, defined as phase difference 
along the longitudinal axis between z = 6 fim and z — fim 
(after unwrapping), plotted as a function of time, as the laser 
is turned up to Uf = 40Er in 200 ms. 



the overall phase profile to vary in time. The second, 
more important, source of the oscillations is the induc- 
tion of transverse and longitudinal collective modes of 
the condensate. As the laser is turned up, the transverse 
confinement increases (see Fig. , and the BEC is pulled 
inward in the transverse direction and pushed outward in 
the longitudinal direction. This flow of supcrfluid gives 
rise to density oscillations along the longitudinal and ra- 
dial axes, as illustrated in Fig.^J As the lattice deepens, 
the frequency of the axial oscillations decreases because 
of the increasing effective mass |60j and the time between 
classical turning points (where the phase is most flat) 
lengthens. 

For a sufficiently deep lattice, the axial dipole mode 
period exceeds experimental timescales and the accumu- 
lated phases become effectively locked in. If the radial 
modes in each well all have the same frequency and are in 
phase, then the axial phase profile would be unaffected by 
their presence. In fact, for non-adiabatic ramps each dis- 
connected well has a slightly different radial confinement 
frequency, so one would expect the longitudinal phase 
profile to change with time even for deep lattices. As dis- 
cussed in the next section, these phase oscillations give 




FIG. 12: Density oscillations, (a) Density (in units of fim~ 3 ) 
along half the longitudinal axis as a function of time as laser 
is turned up to 40Er. (b) Density (in units of /im" 3 ) along 
half the transverse axis as a function of time as laser is turned 
up to 40-Efl in 200 ms. On the whole the BEC is getting wider 
in the longitudinal direction and narrower in the transverse 
direction, but the width oscillates in both directions along the 
way. The oscillations are very large for about the first 100 ms, 
and much smaller after that. 



rise to 'collapses and revivals' of the interference pattern 
which were not observed in the OTFYK experiments. 



VIII. ZETA AND VISIBILITY 

For more quantitative comparisons with experiment, 
it is useful to calculate the quantities used by OTFYK 
to compare their number squeezing model, their effective 
one dimensional GP model, and their data, namely £ and 
visibility jEJ [3 . £ is defined as the ratio of the width 
of a single peak to the distance between the peaks and 
Co = 0.22 is the value of £ for Uf = 6Er. To determine 
C, we fit the cross section of the smoothed absorption 
probability through the longitudinal axis z, A(0, z), to a 
double Gaussian: 

Be -( z - Zl f/2a* +Ce -(z-z 2 f/2*z (?) 
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FIG. 13: Smoothed absorption cross sections (dots) and dou- 
ble Gaussian fits (solid lines) for a range of laser strengths. 
The fits are not perfect, especially at high laser strengths, 
but the differences between the fits and the data are much 
smaller than the differences between the data at different laser 
strengths. The fit shown for 40Er is the worst fit in the data 
set. 

and then £ is the ratio of the width of the Gaussians to 
the distance between their centers: 

C = a/(zi - z 2 ) (8) 

and Co is the value of £ for Uf = 6Er. Fig. 1131 shows a 
sample of smoothed absorption cross sections and double 
Gaussian fits for a range of laser strengths. Visibility 
is defined as the difference between the average of the 
maxima of the two peaks and the minimum between the 
peaks, divided by their sum. 

In a further analysis of the OTFYK experiment, Tuch- 
man [53j plots £ and visibility as a function of Ngfi/j, 



where N is the is number of particles in two wells, g is 
defined as in Eq. (JIJ, and (5 and 7 are defined by the 
following integrals over localized wave functions: 

13 = J d\^(v) (9) 

1 = f d 3 rMr){-^ 2 + V ext (r)}Mr) (10) 
J 2m 

Unlike f and visibility, which are determined by fits 
to experimental data, the parameter Ng/3/j is derived 
from a two-mode model in which there are two fixed wave 
functions, ipi and ip2 , in each of two potential wells in the 
optical lattice. This parameter is essentially a measure of 
laser strength and density, and it is approximately expo- 
nential in laser strength. Our own calculations show that 
this factor varies somewhat depending on the details of 
the theoretical model used to determine it. Therefore, we 
prefer to plot C and visibility versus laser strength Uf, an 
experimentally determined parameter, rather than versus 
Ngfi/j, which is based on particular theoretical model. 
According to the numbers given in Ref. [3^|, a range of 
laser strengths from &Er to bQE^ corresponds approxi- 
mately to a range of NgPh from 10 1 / 2 to 10 5 . 

Fig. 1141 shows C an d visibility for our simulations as a 
function of Uf (laser strength upon release). Each point 
corresponds to a simulation in which the laser strength 
is raised to Uf in 200 ms, the magnetic trap is turned off 
in 40 /is, the BEC is held in the laser for 2.785 ms, the 
laser is turned off and the BEC is allowed to expand for 
8 ms. For each point, £ and visibility are calculated from 
a cross section of the smoothed absorption probability 
discussed in Sec. IIV CI 

Perhaps the most striking aspect of the numerical re- 
sults is that the values of both £ and visibility do not 
merely increase or decrease as a function of Uf, but os- 
cillate. These oscillations are due to the excitation of 
collective modes as discussed in the previous section, and 
do not appear in the experimental data or in the models 
employed in Refs. |38|,|5^]. It is possible that the oscilla- 
tions present in the GPE simulations would be strongly 
damped in actual experiments, with the high radial ener- 
gies being transferred to high-lying axial modes. When 
atoms are loaded into blue-detuned lattices, where the 
transverse confinement due to the external magnetic trap 
is weak compared to that of red-detuned lattices, consid- 
erable radial heating has been observed leading to 
loss of interference contrast 0] . 

The observation of phase oscillations in numerical sim- 
ulations leads to predictions that would be interesting 
to test experimentally. If the BEC were released when 
there was a peak in the amplitude of the phase oscilla- 
tions, the interference pattern could be lost for a rela- 
tively low barrier height. If one waited a little longer to 
release the BEC until the phase flattened out again, the 
interference pattern would return for a few oscillations 
until disappearing due to damping. Note that these 'col- 
lapses and revivals' are completely driven by mean- field 
effects, and are unrelated to similar phenomena found 
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in 3D lattices [64|. Such effects were not reported in 
the OTFYK work, but contrast oscillations attributed to 
quantum fluctuations were recently observed by the same 
group after non-adiabatically rampinga deep lattice to a 
relatively shallow depth of 16.6-Er |54j . 

While we cannot compare the results shown in Fig. 1141 
directly to the experimental results because we are plot- 
ting our results as a function of Uf rather than Ng(3/^f, 
we can make a qualitative comparison based on the given 
ranges of Uf and Ng/3/-f. Our results are qualitatively 
similar to those found in experiments in that the 
fringe contrast decreases approximately linearly with lat- 
tice depth. With our best estimate of the correspondence 
between Ngf3/j and Uf, the average values of ( found in 
the simulations are consistently smaller (on the order of 
30%) than the experimental data. Besides this uncer- 
tainty, this systematic difference is likely due to at least 
two factors. First, the mean-field calculations cannot in- 
corporate the physics of number squeezing that is likely 
to be a contributing factor in real systems. Second, the 



damping of the induced collective modes discussed above 
should diminish the effective scatter of the data points (of 
order 25% of the mean values shown in Fig. lTlfl . heating 
the system and reducing the overall interference visibility. 



IX. CONCLUSIONS 

We have shown that mean-field effects, as modeled by 
the three-dimensional time-dependent GPE, can explain 
both the loss and return of interference for a BEC in 
a one-dimensional optical lattice. The simulations yield 
behavior that is qualitatively similar to that observed of 
the OTFYK experiments [38| without the need to in- 
voke physics beyond mean-field theory, such as number 
squeezing or condensate fragmentation. The central re- 
sult of the present work is that ensuring adiabaticity dur- 
ing the loading of a BEC into a deep optical lattice is 
experimentally difficult to achieve, and that neglecting 
the rather large effects of mean-field excitations will lead 
to an incomplete description of future experiments. 

The results of GPE simulations differ from the exper- 
imental data in a few small respects, however. The loss 
of contrast in the interference patterns, after dropping a 
deep optical lattice and allowing the cloud to expand, is 
not quite as pronounced as in the experiments. One ex- 
planation for the small difference is likely to be the pres- 
ence of number squeezing, which is not captured by the 
GPE. Another possibility is that the radial excitations 
that are induced by the lattice ramp will rapidly damp 
into high-lying axial excitations in real experiments (the 
GPE has no damping mechanism), leading to heating 
and its attendant loss of contrast. These collective modes 
lead to phase oscillations that in turn yield periodic vari- 
ations in the fringe contrast; though these 'collapses and 
revivals' were not seen in the original OTFYK experi- 
ments, similar oscillations attributed to quantum fluctu- 
ations have been observed more recently by the members 
of the same group [54| . Because the creation and destruc- 
tion of interference patterns are widely used techniques 
to measure the presence of coherence and its loss, it is im- 
portant to understand all the mechanisms that can affect 
this important experimental measure. 
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APPENDIX A: NUMERICAL METHODS 

The numerical simulation of the full three dimen- 
sional time-dependent Gross-Pitaevskii Equation (GPE) 
is performed in Cartesian coordinates using a C code 
adapted from earlier work in the Reinhardt group. The 
time integration uses the variable step fourth-fifth order 
Runge-Kutta integrator odeint from Numerical Recipes 
in C (g^. The spatial integration uses a pseudospectral 
method l6al with fast Fourier transforms from the fftw 
library [67]]. The basic idea behind the pseudospectral 
method is that the wave function is expanded in terms of 
coordinate discretized trigonometric functions, reducing 
a partial differential equation into a set of coupled ordi- 
nary differential equations via fast Fourier transforms to 
coefficient basis. We have adapted the code using Type II 
Fourier transforms, which are appropriate to the bound- 
ary conditions of half of a symmetric box, to take ad- 
vantage of the symmetry of the problem and reduce the 
analysis to one quadrant of the system. This adapta- 
tion increases the speed of the calculations by a factor 
of 8. In theory, since the system of interest is cylindri- 
cally symmetric, further improvement in speed could be 
achieved by using cylindrical coordinates and thus effec- 
tively reducing the problem two dimensions. In prac- 



tice, the speed of the fast Fourier transforms appropriate 
to Cartesian coordinates eliminate the disadvantage of 
working in 3D and the 3D Cartesian code runs signifi- 
cantly faster than a 2D cylindrical coordinate code using 
the discrete variable representation method rather than 
the pseudospectral method. 

For the results shown in this paper, we use 24.86 grid 
points per micron in the longitudinal direction (10.44 grid 
points per site) and 1.09 grid points per micron in the 
transverse directions. During the expansion of the BEC, 
we reduce the number of grid points per unit length in 
the longitudinal direction by a factor of 2 during the first 
4 ms, and then by an additional factor of 2 during the 
next 4 ms. We have done sample runs with up to twice 
as many longitudinal grid points and four times as many 
grid points in each transverse direction and checked that 
this does not change the results. Many more grid points 
are needed in the longitudinal direction than in the trans- 
verse direction because of there is much more variation 
due to the laser, which has a wavelength of 0.84/im. 

The simulations were done on a computer with four 
3 GHz Xeon processors and 2 Gb of memory running Red 
Hat Linux. On this machine, using the Intel C compiler, 
the computations of the BEC expansion take about 4 
hours using the minimum necessary number of grid points 
and the 200 ms laser turn-on computations take about 
40-60 hours. (The code runs about 1.5 times faster when 
compiled with the Intel C compiler than when compiled 
with gcc.) 
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